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I. Introduction 



Simulations including two degenerate light flavours and a non-degenerate doublet of quarks 
are currently being performed by the European Twisted Mass (ETM) Collaboration. The inclusion 
of rif = 2+1 + 1 flavours is a necessary step to move towards a realistic situation. Fermions are 
described by the maximally twisted mass lattice QCD (MtmLQCD) action [1] and gluons by the 
Iwasaki action [2]. While the first physical results are very encouraging [3], dedicated simulations 
are required to perform the non-perturbative renormalization of operators in a mass-independent 
scheme, where renormalization constants (RCs) are defined at zero quark mass. In the study of rif = 
2+1 + 1 QCD ETMC is adopting the RI-MOM scheme [4]. The RCs are evaluated by extrapolating 
to the chiral limit the RC estimators computed in the theory with rif = 4 mass degenerate quarks 
for a range of mass values l . Here we report on the progress we made in this project. 



1.1 Action and quark mass parameters 

For the present study we consider the lattice action 

Sl = Sj™ + a 4y £ t Xf [/• V - |V*V + m , / + /75 wl */(*) d-D 

x f=l 

where is a one-flavour quark field in the so-called twisted basis and in this work rf is set to 
either 1 or —1. Passing from the twisted to the physical quark basis 2 

S L = S™ + a 4 £ £ q f \y- V - ^^./(-f V*V + m cr ) +M 0J \ g f (x) . (1.3) 

x f=i 

The bare mass parameters can be rewritten as 



I Mo f — YYIq^ LI /• 

Mq j = yj {mo j - m cr ) 2 + fij , sin G 0J = ^— cr , cos G j = . (1.4) 



Their renormalized counterparts read M f = Z P M f = yJz\ m lcAC + M/ and tan©/ = Za ^ cac . The 
parametrization in terms of M and 9 is convenient because the leading term of the Symanzik local 
effective Lagrangian involves only M, not G. As we will see later (see the end of section 2), this 
remark is at the basis of our method to obtain 0(a) -improved RC-estimators at all scales even out 
of maximal twist. Since, for practical reasons, we work in a partially quenched setup with all four 
flavours having equal mass parameters, we will have to consider in our analysis four quark mass 
parameters: M sea , se a, A^vai, Gvai- 



^or Monte Carlo simulations we used a highly optimized implementation of a HMC-like algorithm [5]. 
2 The relation between twisted (Xf fields) and physical (qf fields) quark basis is 

X f ^q f = exp[|(f - Ooj)y 5 r f }Xf, Xf -> Qf = Z/exp[|(f - f)y 5 r f } (1.2) 
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1.2 RF-MOM scheme and our setup 



The focus of the present study is on flavour non-singlet quark bilinear operators, Or = Xf^Xf 
(or %fiT%f), with T = V,A, T, which are written in terms of X and X quark fields (i.e. in the 
standard quark basis for untwisted Wilson fermions). RCs are named after the expression of the 
operators in this basis so as to match the usual notation in the literature about Wilson fermions. 
As convenient in lattice studies, we adopt the RF-MOM scheme [6, 4], which is defined as follows. 
A first condition fixes the quark field renormalization, namely 



Z" 1 



\2N{p) L 



Mr P s f (p)- 1 ) 
Pp 



1, any/, 



(1.5) 



where p 2 = , Pn = ^ sin ap^. The sum Yip only runs over the Lorentz indices for which p p 

is different from zero and N(p) = Ep 1 • The renormalization condition for the operators Or reads 



A^ ff '\p,p)Pr 



= 1, S + f- (1-6) 

Above Sf(p) — a 4 Y* x e~ ipx (Xf( x )Xf(ty) 1S the Xf Held propagator in momentum space, while 

4 f \p,p) = Sj l (p)G { r ff '\p lP )S-/(p) (1-7) 
denotes the quark bilinear vertex that is obtained by "amputating" the Green function 

G ( / f '\p,p) = a 8 £ e -<>(^) (xf(x)(xfrx f >)(o)Xf>(y)) r = s,p,v,a,t. U-8) 



Barring cutoff effects, RCs are independent of sign(rf). For practical reasons here we limit our- 
selves to rf = —rf in evaluating Zp, see eq. (1.6). 



2. Strategy for RCs in the N f = 4 theory 

In order to extract useful information from simulations performed with twisted mass Wilson 
fermions one must know the twist angle, CO = | — 0, with good precision. The level of precision 
requested for CO depends on the observable of interest. In our case, after an exploratory study on a 
few 16 3 x 32 lattices [7], and some tests near maximal twist on a 24 3 x 48 lattice we have chosen 
to work out of maximal twist. 

Figure 2 illustrates the difficulties of tuning to maximal twist, i.e. setting mpcAC to zero, in the 
simulation setup for RC computations, at least if the lattice spacing is not very fine. Specifically, 
the slope of mpcAC vs 1/ (2k) in figure 1(a) suggests that near mpcAC = simulations are in a region 
with a sharpe change of the slope for mpcAC where it is difficult to extract useful information. On 
the other hand figure 1(b) gives a more quantitative view of this problem showing results from one 
simulation close to the critical point (the point closest to mpcAC = in figure 1(a)). It appears that 
due to the long fluctuations a precise measurement of the PCAC mass will require for this case a 
very large number of Monte Carlo trajectories. In fact, we have observed a similar feature for all 
the ensembles with \am FC Ac\ < 0.01 at both /3 = 1.95 and j8 = 1.90 (i.e. for a > 0.08 fm). 
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In summary, working at maximal twist for the chosen range of twisted masses (see table 1) 
would imply a considerable fine tuning work owing to the difficulties in determining am?cAC near 
ampcAC = 0. To alleviate the problem one would need to increase the value of the twisted mass, 
jif, and thus Mf. Instead, working away from maximal twist, one can avoid the metastable region 
of parameter space and measure the twist angle with good precision. This comes at the price of 
a moderate increase of the quark mass Mf and of a slightly more involved analysis. In our RC- 
estimators cutoff effects linear in a are expected to be small and can anyway be removed with 
controlled precision by averaging the results obtained for a given Mf at opposite values of Of. 



&< 0.00 




3.090 3.095 3.100 3.105 3.110 3.115 

1/2k 




(a) 



01 fo00 1500 2000 2500 3000 3500 4000 

MC time 
(b) 



Figure 1: (a) am?cAC versus 1/2k on lattices 24 3 x 48 at j3 = 1.95; (b) Monte Carlo history of am?cAC f° r 
the most critical case (corresponding to the red cross point) in panel (a). 

In fact, from the symmetry of the lattice action Sl under x (0 O — >> — 0o) x^x (Mo — >> 
—Mo) [1, 8, 12] it follows that the 0(a 2k+l ) artifacts occurring in the vacuum expectation values 
of (multi)local operators O that are invariant under & x (0q — > — Oo) are quantities that change 
sign upon changing the sign of 0q (or 0). Hence 0{a 2k+l ) cutoff effects vanish in 0-averages: 



(0)W e + (0)\ 



M,-G 



The same is true for operator form factors invariant under x (0q -> — %) 



and, in particular, for our RC-estimators at all values of Mf and p 2 . 



3. Current analysis and preliminary results 

Here we detail the analysis procedure we followed in order to obtain very preliminary results 
on the RCs of interest. Indeed, at this stage our main goal was checking the feasibility of the 
project. In particular, the analysis procedure is not yet the optimal one, for instance concerning 
the order of the various steps, and some refinements, such as the subtraction of the known cutoff 
effects at 0(a 2 g 2 ) [9], are still omitted. While these improvements will be included in the final 
analysis, the present work shows that the strategy advocated in section 2 allows to extract the RCs 
of the quark field and quark bilinear operators with a ~ 1% level precision by means of stable 
simulations at a lattice spacing (a ~ 0.08 fm) which is among the coarsest ones explored in the 
study of rif = 2 + 1 + 1 QCD by ETMC. 

In practice, for a sequence of M sea -values, we produced for each M sea two ensembles with 
opposite values of sea . We label them as Ep/m, where E = 1,2... and p/m refers to sign(0 sea ). On 
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ensemble 


^Msea 






gsea 


ajl 


val 




lm 


0.0085 


-0.04125(13) 


0.03288(10) 


-1.3093(8) 


[0.0085,. . 


. , 0.0298] 


-0.0216(2) 


lp 


0.0085 


+0.04249(13) 


0.03380(10) 


1.3166(7) 


[0.0085,. . 


, 0.0298] 


+0.01947(19) 


3m 


0.0180 


-0.0160(2) 


0.02182(9) 


-0.601(6) 


[0.0060,. . 


, 0.0298] 


-0.0160(2) 


3p 


0.0180 


+0.0163(2) 


0.02195(9) 


0.610(6) 


[0.0060,. . 


, 0.0298] 


+0.0162(2) 


2m 


0.0085 


-0.02091(16) 


0.01821(11) 


-1.085(3) 


[0.0085,. . 


. , 0.0298] 


-0.0213(2) 


2p 


0.0085 


+0.0191(2) 


0.01696(16) 


1.046(6) 


[0.0085,. . 


. , 0.0298] 


+0.01909(18) 


4m 


0.0085 


-0.01459(13) 


0.01409(8) 


-0.923(4) 


[0.0060,. . 


, 0.0298] 


-0.01459(13) 


4p 


0.0085 


+0.0151(2) 


0.01441(14) 


0.940(7) 


[0.0060,. . 


. , 0.0298] 


+0.0151(2) 



Table 1: Mass parameters of the ensembles analysed for this contribution. From the formulae in sect. 1.1 it 
follows that in the valence sector we have 0.013 < aM val < 0.033 and 0.4 < 1 va i | < 1 .2 (Gvai/wpcAC > 0). 




'••'•V •■ 

(a) (b) 



Figure 2: For the example of ensemble 4m, (ap) 2 — 1.5: (a) subtraction of Goldstone pole contribution and 
valence chiral extrapolation in Tp = Tr [ApPp]; (b) overview of the valence chiral extrapolation for all RCs. 

each ensemble Ep/m, with (M^ m , 0s&J m ) we compute the RC-estimators for several values of the 
valence mass parameters (M va i, va i) and p 2 (all corresponding to "democratic" momenta p, in the 
sense specified in [10]), as summarized in table 1. Then we proceed in various steps as follows. 

Valence chiral limit. A fit of RC-estimators linear in (M^j) 2 turns out to be numerically adequate 
(see fig. 2(b)). For r = P (see fig. 2(a)) or, due to 0(a 2 ) terms, r = S, we have also kept into 
account the contribution oc (M^})~ 2 coming from the Goldstone boson pole. 

0(a 2 p 2 ) discretization errors. We applied two different methods, following [10]. In the first 
method ("Ml"), after bringing, via the known [11] perturbative evolution the RC-estimators 
to a common renormalization scale {py^ = I /a 2 ), we remove the remaining 0(a 2 p 2 ) dis- 
cretization errors by a linear fit in p 2 . The second method ("M2") consists in simply taking 
the value of the RCs estimators at a high momentum point fixed in physical units. We chose 
= 12.2 GeV 2 . The two approaches yield RC results differing only by cutoff effects. 
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Figure 3: Residual /? 2 -dependence of RC-estimators at scale 1 /a and RC values from method Ml for the 
cases of ensemble lm (panel (a)) and 2p (panel (b)). 
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Z q (Ml) 






Z P (M1) 
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* 


Z 4 (M1)-0.1 


1 


♦ 


Z S (M1) 


♦ 


♦ 


Zy(VJ\)-Q.l 



0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 

(aM) 2 



(a) 



<3 0.6 



♦ ♦ Zq(M2) 

♦ ♦ Za(M2) 

Zp(M2) 

♦ ♦ Zs(M2) 



0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 

(aM) 2 



(b) 



Figure 4: M sea -dependence before (empty symbols) and after (full symbols) 6 -average of RC-estimators 
for a few operators as well as their chiral limit value. The left (right) panel corresponds to Ml (M2) results. 
In (a) we also show Zy(WI), which is obtained by exploiting an exact lattice Ward-Takahashi identity (WI). 



Removal of 0(a) artifacts. It is achieved by 0-average (see section 2) of the RCs estimators, 

1 



Zr(M s E ea ,|0 } 



,~ E I 

sea' I ^seal 



(3.1) 



where parameterizes the dominating 0(a) effects in RC-estimators that (in the present 

analysis) arise from employing A^. Ep ( m ) in the valence chiral extrapolation. 

Sea chiral limit. The quantities Zr (M^ ea , 1 S g a | ) are extrapolated to M sea = by using the fit Ansatz 



Zr(Msea,0s< 



Z r + AM s 2 ea + £M s 2 ea cos(2ft 



(3.2) 



This Ansatz can be justified by an analysis a la Symanzik of the lattice artifacts in Zr(M sea , 9 S{ 
up to 0(M 2 ea ) and neglecting chiral spontaneous symmetry breaking effects [12]. 



The first, very preliminary results of this analysis are summarized in table 2. 



6 



RCs - four dynamical flavours 



David Palao 



Method Z A Z v Z P (l/a) Z s (l/a) Z P /Z S Z T (l/a) Z q (l/a) 



Ml 0.761(08) 0.630(05) 0.438(08) 0.614(09) 0.716(21) 0.753(07) 0.767(06) 
M2 0.771(03) 0.674(03) 0.496(04) 0.647(03) 0.767(08) 0.768(03) 0.813(02) 



Table 2: Preliminary RC results at j3 = 1.95 from the analysis of section 3. We also get Zy(WI) = 0.612(1). 

4. Conclusions and outlook 

We have described our strategy to compute 0(a) improved operator RCs for the Nf = 4 lattice 
action currently used by ETMC. We have shown that the method advocated in this work provides 
very encouraging results at one lattice spacing (a ~ 0.08 fm) that is among the coarsest simulated in 
the study of QCD with rtf = 2+1 + 1 dynamical flavours. In particular, the observed dependences 
of RCs on valence and sea quark masses are mild and quite in line with our experience [10] in 
rif = 2 QCD. Besides the technical improvements mentioned in section 3, we plan to possibly add 
few more ensembles at a ~ 0.08 fm (j8 = 1.95) and to extend our work to other lattice spacings. 

We thank IDRIS and INFN/apeNEXT for giving us CPU time necessary for this study. 
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